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ABSTRACT 

The  long  time  behavior  of  the  solutions  of  nonlinear  parabolic  initial 
value  problems  on  the  line  has  been  investigated  by  many  authors.  In 
particular  they  have  shown,  under  certain  assumptions,  the  existence  of 
traveling  waves  to  which  a  large  class  of  initial  data  eventually  evolves. 

They  have  also  proved  that  which  traveling  wave  solution  is  picked  out  as  the 
asymptotic  state  often  depends  on  the  behavior  of  the  initial  data  at 
infinity.  This  causes  difficulties  for  the  numerical  simulation  of  the  long 
time  evolution  of  such  problems.  In  particular,  if  an  artificial  boundary  is 

introduced,  the  boundary  condition  imposed  there  must  depend  on  the  initial 

- 7  .  . 

data  in  the  discarded  region/  In  this  work  we  derive'- such  boundary 

—  r 

conditions,  based  on  the  Laplace  transform  solution  of  the  linearized  problems 

r  S 

at  -Wtf  illustrate  their  utility  by  presenting  a  numerical  solution  of 

Fisher's  equation,  a  nonlinear  parabolic  equation  with  traveling  wave 
solutions  which  has  been  proposed  as  a  model  in  genetics. 
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SIGNIFICANCE  AND  EXPLANATION 


Nonlinear  partial  differential  equations  of  parabolic  type  arise  in 
various  applications.  Examples  include  models  of  chemical  kinetics  and 
population  dynamics.  In  many  cases  it  is  the  evolution  of  the  initial  data 
into  some  simple  final  state  which  is  of  interest.  For  a  class  of  initial 
value  problems  on  the  line,  other  authors  have  shown  that  the  final  state  is 
usually  a  traveling  wave  and  is  determined  by  the  initial  data  at  infinity. 
In  this  work  we  present  a  method  for  the  numerical  simulation  of  this 
evolution. 


As  a  finite  domain  is  required  for  the  numerical  method,  it  is  necessary 
to  introduce  artificial  boundaries.  The  boundary  conditions  imposed  there 
must  depend  on  the  initial  data  in  the  discarded  regions  if  the  correct  long 
time  solution  is  to  be  found.  We  construct  such  conditions  using  the  Laplace 
transform  solution  of  the  linearized  problems  at  ±“.  Their  utility  is 
illustrated  by  the  solution  of  Fisher's  equation,  a  model  of  the  spatial 
advance  of  an  advantageous  gene.  It  is  hoped  that  this  method  will  give 


THE  NUMERICAL  CALCULATION  OF  TRAVELING  WAVE  SOLUTIONS 
OF  NONLINEAR  PARABOLIC  EQUATIONS  ON  THE  LINE 


1  INTRODUCTION 


Thomas  Hagstrom  and  K.  B.  Keller* 


We  consider  the  numerical  solution  of  the  Cauchy  problem  for  a  class  of  nonlinear 
parabolic  equations* 


a) 

“t  *  f(uxx'ux'u) 

.  —  <  x  <  «, 

t  >  0  ; 

b) 

lim  u(x,t)  “  ^+, 

X-«» 

lim  u(x,t)  “ 

1 

c) 

u(x,0)  -  Uq ( x )  > 

d) 

3 

3^  f(a,b,c)  >  1 

for  all  a,b,c 

. 

We  assume  that  satisfy: 


(1.2) 


f<0,0,$±)  -  0  I 


fu<°,°>±)  *  0  f 


and  that  the  initial  data.  u0(x),  satisfies  (1.1b).  In  particular,  we  are  Interested  in 
simulating  the  evolution  of  the  initial  data  into  traveling  waves. 

Hagan  [3,4]  has  presented  an  extensive  analysis  of  problem  (1.1).  We  paraphrase  some 
of  his  results  below: 

(i)  Nonmonotonic  waves  are  unstable  in  general. 

(ii)  The  stability  of  monotonic  waves  of  speed  c  can,  in  general,  be  determined  by 
an  examination  of  their  trajectories  in  the  phase  plane  of: 

W  ’“VI 

(1.3) 

f(v',v,w)  +  cv  -  0  . 

(ill)  If  traveling  waves  exist,  a  large  class  of  initial  data  satisfying  (1.1b)  will 
evolve  to  a  traveling  wave. 


(iv)  In  certain  situations,  infinitely  many  wavespeeds,  c,  are  allowed.  In  this 
case,  the  traveling  wave  which  is  eventually  seen  depends  on  the  behavior  of  the  initial 
data  at  infinity. 

The  numerical  solution  of  (1.1)  requires  a  finite  computational  domain.  One  way  to 
obtain  such  a  domain  is  to  introduce  artificial  boundaries  at  the  points 
x  “  r+  >  t_,  and  impose  boundary  conditions  there.  The  derivation  of  proper 

boundary  conditions  is  the  main  subject  of  this  work.  A  general  theory  of  boundary 
condition,  at  an  artificial  boundary  is  given  by  the  authors  in  [6].  This  theory  is  not 
directly  applicable  to  time  dependent  problems  in  unbounded  spatial  domians  such  as 
(1.1).  However,  a  Laplace  transformation  in  time  yields  a  problem  of  the  right  form.  In 
section  2,  conditions  are  derived  for  the  transformed  problem  and  inverted  to  yield 
conditions  in  the  real  variables. 

We  note  that  use  of  the  proper  boundary  condition*.  is  crucial  whenever  (iv)  holds. 
Then,  the  "naive"  conditions: 

u(T+,t)  -  *+  , 

(1.4) 

U(T_,t)  -  , 

must,  in  general,  fail  to  lead  to  the  correct  long  time  solution. 

In  section  3  a  specific  problem  of  the  form  (1.1)  is  introduced:  the  Cauchy  problem 
for  Fisher's  equation.  It  has  traveling  wave  solutions  of  all  speeds  c  >  2.  Gazdag  and 
Canosa  [1]  present  a  numerical  solution  of  Fisher's  equation  using  boundary  conditions 
analogous  to  (1.4).  As  predicted  by  the  theory,  their  solution  always  evolved  to  the 
traveling  wave  of  minimum  speed.  Here  we  present  calculations  using  the  boundary 
conditions  derived  in  section  2.  The  numerical  solution  is  seen  to  evolve  to  the  correct 
traveling  wave  for  a  variety  of  choices  of  initial  data. 

We  note  that  the  method  of  deriving  boundary  conditions  presented  here  can  be  applied 
to  other  time  dependent  problems,  including  some  problems  of  hyperbolic  type.  For  other 
examples  the  reader  is  referred  to  Gustafsson  and  Kreiss  [2]  and  Kagstrom  [5]. 

The  authors  thank  Prof.  J.  D;  Murray  for  bringing  this  problem  to  their  attention. 


2.  CONSTRUCTION  OF  THE  BOUNDARY  CONDITIONS 

We  now  construct  boundary  conditions  at  the  right  boundary,  x  *  t+.  (The 
construction  at  the  left  will  be  analogous-)  Only  the  linearized  problem  in  the  tail  is 
considered  and  a  coordinate  system  moving  to  the  right  with  speed  c  is  assumed; 

a)  v(x,t)  “  u(x,t)  -  x  >  t+  ; 

b)  vt  "  f1vxx  +  f2vx  +  cvx  +  f3v  * 

(2.1) 

c)  v(x,0)  **  Ug(x)  -  $+  > 

d)  lim  v(x,t)  =  0  ; 
x-*“ 

where  the  constants  f ^  are  given  by: 

fi  -  irs-r  (0'°'V  ' 

XX 

(2.2,  f2  “  ff=TT  (0'0'V  ' 

X 

f3  -  !£  (0,0, *+)  . 

Following  the  general  ideas  presented  by  the  authors  in  [6] ,  two  problems  must  be 
solved;  boundary  conditions  for  the  homogeneous  problem,  (2.1b,d)  combined  with  zero 
initial  data,  must  be  found  as  well  as  a  particular  solution  of  (2.1b,d)  which  satisfies 
(2.1c).  The  homogeneous  problem  is  considered  first. 

Boundary  Conditions  for  the  Homogeneous  Problem 
We  introduce  the  temporal  Laplace  transform: 


w(x,s)  “  /  e  ®  ui(x,t)dt 
0 


If  <i>  is  a  solution  of  (2.1b,d)  with  zero  initial  data,  then  w  satisfies: 


f„u>  +  (f_  +  c)»  +  (f,  ~  “  0  t 

1  xx  2  x  3 


lim  u(x,s)  “  0 


Equation  (2.3a)  has  the  basic  exponential  solutions: 
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*  k  (s;c)x 
01  -  e  1 


where  k^(sic)  are  given  by: 


~^2  +  1  2  1/2 

kt  ■  — 1 2?;  (<f2  ♦  c)  +  4f,<8  -  V1 


For  real  part  of  a  sufficiently  large,  k+  will  have  real  part  positive  and  k_  will 
have  real  part  negative.  (Recall,  f^  >  1.)  Hence,  the  admissible  solution  has  exponent 


k_.  It  satisfies : 


which  can  be  rewritten* 


U>x(T  +  |S)  “  k_(8fC)ttl(T  +  |s)  ! 


U  <T  +  JS)  a 

k'Ts,cl  *  “(Va>  * 


Using  the  convolution  formulas  and  the  expression  for  the  inverse  transform  of  —  (see, 
e.g.,  Oberhettinger  and  Badii  [8]),  (2.7)  can  be  expressed  in  the  real  variables: 

_  t  tf  -a2](t-p)  2  _ 

-/f.  /  e  f  -  -  «ea  p  Erfc(a/t  -  p)]<o  (t  ,p)dp  -  id(t  ,t)  ; 

0  A(t  -  p)  x 

(2.8) 

fj  ♦  c 
a  ;  -  . 


Particular  Solution 

We  now  find  particular  solutions,  assuming  that  u0  -  $+  can  be  expressed  as  a  finite 
sum  of  exponentials: 


N  -a  <x-x+) 

u0(x)  -  *+  -  I  a  e  .  a  >  0,  x  >  t+  . 

■*»1  J  ■> 


We  note  from  Hagan's  analysis  [3]  it  is  necessary  in  many  cases  that  uQ  -  $  decay  at 
least  exponentially  if  traveling  wave  solutions  are  to  be  found.  From  (2.5),  with 


-4- 


8  -  0,  we  see  that  traveling  waves  of  speed  c  have  exponential  decay  rates  given  by: 


From  (2.9)  we  see  that  given  any  exponent,  -cr,  there  exists  a  unique  speed  such 


(2.10) 


It  is  given  by: 


-“j  =  9_(Cj)  or  -a..  =  g+Jcj)  . 


(2.11) 


'  “i  +  '  f2 


Hence,  each  exponential  can  he  associated  with  a  unique  traveling  wave  solution,  from 
which  a  particular  solution  can  be  found: 


(2.12) 


N  a.T  -a  (x-(c.-c)t) 

V  *  -  3  J  3 


V  (x,t)  -  l  a  e  3  e 

P  j-1  3 


Combining  (2.8)  and  (2.12)  yields  the  following  linearized  boundary  condition  for  u 


at  x  ■  x+: 


(2.13) 


__  t  [f  -o‘](t-p)  ,  2.  .  _  N 

-/f.  /  e  f  -  ~  “e  p  Erfc(aZt-p)  ]  (u  (t  ,p)  +  £  a  d 

0  /x(t-p)  j-1  3  3 

N  a. (c.-c)t 

-  u(T+,t)  -  *+  -  l  d  e  3  t 


N  a  (c.-c)p 

l  3  3  )dP 


f2  +  c 

a)  a  -  - - : 


(2.14) 


N  -a  (x-x+) 

b)  Uq ( x )  -  ♦.  +  \  d  e 

j-1  3 


*  *  > 


Conditions  at  the  Left  Boundary 

A  similar  boundary  condition  can  be  derived  at  the  left  boundary,  x=  r_.  in 
transform  variables,  a  solution  to  the  linearized,  homogeneous  problem  on  (-“,r  ]  must 
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satisfy s 
(2.15) 


<0  (T  ;  s ) 
x 

k+(s»c) 


U»(T  (S)  , 


(2.16) 


”^2  +  1  2  -  1/2 

- - -  +  ~  Uf,  +  c ) *  +  4f(s  -  f_»  1/2 

2ft  2r1  2  13 


He  now  have: 


fi  “  17T7  (0'°'*-)  ’ 

XX 


(2.17) 


f2  "  57TT  (0'°'*-)  1 

X 


?3  mU  <°'°'+-) 


The  inverse  transform  of  (2.15)  is  given  by: 


/7“  t  (?--«  J(t-p)  1  ^  —2  _  _ 

/  f1  /  e  [  -  —  +  Se  p,Erfc(-a/t  -  p)]w  (T  ,p)dp  -  u»(t  ,t)  ; 

0  /if  (t  -  p)  x 


(2.18) 


f2  ♦  c 


To  find  a  particular  solution  we  assume  that 


M  Cl  (x-T  ) 

V*>  -  ^  5je  ~  '  “j  >  °'  *  <  T. 

j-1  3  3 


Each  exponent,  can  be  uniquely  associated  with  a  linear  traveling  wave  of  speed  c 


through  equation  (2.16)  (with  s  “  0): 


(2. 19) 


C,  *  -f,  -  f.cr  -  —  , 


leading  to  the  particular  solution: 


(2.20) 


C/t )  -  ♦  +  l  d  < 

1-1  3 


M  -a  t  a  (x+(c-c.)t) 
y  d.e  j  %  j  j 


Combining  (2.20)  with  (2.18)  yields  a  linear  boundary  condition  at  T_,  analogous  to 
(2.13); 


—  t  [f  -a2] (t-p) 

Me 


t  -  a2ft-  t _  M  a  (c-c  ,)p 

f -  +  ae  p,Erfc(-a/t-p)Ku  (t  ,p)  -  )  a  .d^e  3  3  )dp 

/w(t-p)  X  j-1  3  3 


(2.21) 


M  _  a  (c-c.)t 
u(t  ,t)  -  ♦  -  J  d.e  3  3 

j-1  3 


-  f2  +  C 
a)  a  -  -  ; 


(2.22) 


M  _  a  (x-t_) 

u0(x)  =  $_  +  I  d  e  ,  x  <  t_ 

1-1  J 


3.  Application  to  Fisher's  Equation 

We  now  apply  the  results  of  preceding  section  to  Fisher's  equation: 


a) 

ut  “  uxx  +  U(1  - 

u),  x  e  (-«• 

(3.1) 

b) 

lira  u(x,t)  -  0, 

lim  u(x,t)  ■ 

x+“ 

x+-» 

c) 

u(x,0)  =  U0(x)  . 

Problem  (3.1)  has  arisen  as  a  model  of  the  propagation  of  an  advantageous  gene.  For  a 
discussion  of  this  application  see,  for  example,  Moran  [7].  it  is  a  special  case  of  (1.1) 
and  various  statements  concerning  the  behavior  of  its  solution  are  consequences  of  Hagan's 
[3]  general  analysis: 

(i)  There  exist  traveling  wave  solutions  of  all  wavespeeds  c  >  2. 


(ii)  All  positive  Initial  data,  Uq(x),  decaying  at  least  exponentially  as  x 


i.*  •. 

i,'  \ 


v. 


.V- 


r:-- 


r.  - 


r.  * 

■\”w 


•  ^ 


VV 


tw* 


& 


■ .  • 


■V 

i 


v:- 


,  •» 
** . 

5? 


v. 

# 

v> 


evolves  to  a  unique  traveling  wave. 


-$x 

(ill)  If  Ug(x)  ~  e  as  x  ♦  •  then  the  solution  evolves  to  a  wave  of  speed 


c(8)  given  by: 


(3.2) 


c(8 ) 


1  +  B 


3 

2 


8  <  1 

6  >  1  . 


The  linearized  boundary  conditions,  (2.13)  and  (2.21),  are  easily  specialized  to  this 
problem.  As  in  section  2,  we  introduce  a  coordinate  system  moving  to  the  right  with 
speed  c  and  choose  and  x_  as  our  artificial  boundary  locations,  we  assume  ♦ 

Uq ( x )  can  be  represented  as  a  finite  sum  of  exponentials  in  the  tails: 

N 


.  -a  (x-t+) 

u0(*>  “  l  die  ,  x  >  r  ; 

1-1  3 


(3.3) 


u.(x)  -  l  d  d 


M  a.(x-T  ) 


+  1,  X  <  T 


1-1 


The  boundary  conditions  we  impose  are: 


a)  -/  e 
0 


t  (i-f-)<t-p>f  ,  c  f(t-p) 


[- 


/ w ( t-p ) 


_  N  ( 1  +a  .  -a  .  c )  p 

Erfc(|  /t-p)](u  (x  ,p)  +  lade  3  3  )dp 

j=1  3  3 


N  ( 1+a^-a  ,c)t 
u(t+,t)  -  l  d^e  33  > 


1-1 


(3.4) 


t  ~(1+  3—) Ct-p)  .  f-<t-p)  _  _  M  (a2+ci.c-1)p 

b)  J  e  f  1  ■  +  5-  e  Erf c( -  r/t-p)](u  (T_,p)  -  j  a.d.e  3  3  )dp 


/ x ( t-p ) 


1=1 


J  J 


(a^tix-l  )t 


u ( t  ,t)  -  1  -  )  d  e  3  3 

1-1  3 
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4  .  1  *  *  m  *  »  ^4  -  s*  •  **  **e*.4.*.*,4.’,-  .  «•  .  •  .  •*,  .*  ,  -*  -  .*  .*  ‘  ■  •  *  -  '  «  4  »  4  . 

»,4V4  •‘-V4  ,v.,'4".^y»  *’•  ,*»  V4  *  •  *  ••  .  .*  •.  V  •  .  . -V-’. 


We  note  that,  by  (3.3),  the  true  solution  should  evolve  to  a  wave  of  speed  c(6)  given  by 


(3.2)  where 

(3.5)  8=  min  {a.}. 

■i  ^ 


In  certain  circumstances  the  particular  terms  in  (3.4a)  have  a  large  exponential  growth  in 
time.  As  this  could  be  a  source  of  error  in  a  numerical  computation,  the  integrals 
involving  them  were  done  exactly.  This  allows  us  to  rewrite  the  right  boundary  condition: 


2  2 
t  f  1-  —  )(t-p)  -j-(t-p) 

!  e  f .  -  f  e  Erfc(|  ✓t-p)]u  (T  ,p)dp 

3  /n(t-p)  X 


(3.4a*  ) 


u(x+,t)  +  l  f  <t>  , 
3  =  1 2  3 


where 


2 

k.  =  1  +  a,  -  a.c  . 

3  3  3 


We  note  that  (3.4a')  explicitly  contains  the  different  evolution  of  initial  data  with  large 
and  small  decay  rates. 

Presented  below  are  the  results  of  some  numerical  computations  of  solutions  of  (3.1) 
using  the  boundary  conditions  (3.4a',b).  A  uniform  grid  was  introduced  and  spatial 
derivatives  were  replaced  by  centered  finite  differences.  The  method  was  implicit  in  time 
and  stable  for  the  ratio  of  the  time  step  to  the  grid  size  sufficiently  small.  At  each 
step  a  nonlinear  system  of  equations  was  solved  using  Newton's  method  with  an  explicit  step 


taken  to  ganarata  tha  initial  guess.  Tha  boundariaa  were  located  midway  batwaan  grldpointa 
and  tha  intagrala  thara  wara  approx imatad  by  tha  traps sold  rula  (away  from  tha 
aingularity) .  for  all  caaaa  daacribad  balow  tha  grid  rangad  batwaan  -12  and  12  and 
contained  171  points.  Tha  time  step  ia  .025,  wall  within  tha  stable  region  in  all  caaaa. 
Initial  conditions  were  generated  in  tha  following  wayt  expanalona  in  tha  tail,  (3.3), 
were  input  and  smoothly  connected  (two  continuous  derivatives)  by  a  combination  of 
polynomial  and  exponential  functions.  The  computations  shown  were  performed  on  a 
VAX  11/780  at  the  University  of  Wisconsin  at  Madison,  though  others  were  done  on  the 
IBM  4341  of  the  Applied  Mathematics  Department  at  the  California  Institute  of  Technology. 

Figure  1  shows  the  evolution,  in  a  coordinate  system  moving  with  speed  4,  of  initial 
data  which  decays,  at  both  ±",  at  a  rate  compatible  with  a  wave  of  speed  4.  The  initial 

data  and  solutions  at  intervals  of  25  time  steps  are  displayed.  A  steady  state  is  reached 

which  must  be  moving  with  speed  4.  Figure  2  contains  the  final  state  (solid  line)  of 

Figure  1.  This  is  the  solution  at  t  »  6.875.  Plotted  with  it  is  a  wave  of  speed  4  found 

by  a  finite  difference  solution  of  the  relevant  steady  problem.  The  agreement  between  the 
two  solutions  is  seen  to  be  excellent. 
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Figure  2 

We  note  that  the  boundary  condition, 

(3.7)  u(x+,t)  »  constant  , 

leads  to  good  results  when  the  speed  of  the  coordinate  system  is  the  same  as  the  speed  of 
the  final  state.  For  a  more  complicated  problem,  however,  this  might  not  be  known  in 
advance.  Indeed,  it  might  the  goal  of  the  computation  to  determine  it.  As  shown  in  Figure 
3,  our  conditions  avoid  this  difficulty.  This  is  the  computed  evolution  in  a  coordiante 
system  moving  with  speed  3  of  the  same  initial  data  used  to  generate  Figure  1.  The  wave  is 
seen  to  move  to  the  right  and,  in  fact,  moves  with  relative  speed  1.  This  is  confirmed  in 
Figure  4,  a  comparison  of  the  solution  at  t  =  6.875  (solid  line)  and  the  wave  of  speed  4 
of  Figure  2  translated  to  the  right  a  distance  of  6.875.  We  believe  the  small  error  at  the 
right  boundary  is  due  to  the  use  of  linearized  boundary  conditions. 

Figure  5  displays  the  computed  evolution  of  initial  data  with  two  decay  rates  in  the 
right  tail:  one  compatible  with  a  wave  of  speed  4, the  other  compatible  with  a  wave  of 
speed  3.  Here,  the  speed  4  part  decayed  at  the  large  rate  while  the  speed  3  part  decayed 
at  the  slow  rate.  As  predicted  by  the  theory,  a  speed  3  wave  is  eventually  reached. 
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Ha  note  that,  a*  It  la  tha  Initial  data  In  tha  right  tall  which  determines  the 
wave speed,  it  la  the  right  boundary  condition  which  la  Important.  Various  choices  for  the 
left  boundary  condition,  for  example  u  “  constant  and  u„  "  0,  ware  tried  and  led  to 
good  results. 

In  susssary,  we  have  shown  that  our  boundary  conditions  consistently  lead  to  correct 
long  tisM  results  while  other  simpler  conditions  do  not.  He  hope  that  their  generalization 
to  sore  complicated  problems,  where  the  final  state  Is  not  known  a  priori,  will  also  give 
reliable  results.  It  should  be  noted,  however,  that  this  has  not  been  proved  even  in  the 
simple  case  described  here. 
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